Scrub encroachment promotes biodiversity in temperate European wetlands under eutrophic conditions

Abstract Wetlands are important habitats, often threatened by drainage, eutrophication, and suppression of grazing. In many countries, considerable resources are spent combatting scrub encroachment. Here, we hypothesize that encroachment may benefit biodiversity—especially under eutrophic conditions where asymmetric competition among plants compromises conservation targets. We studied the effects of scrub cover, nutrient levels, and soil moisture on the richness of vascular plants, bryophytes, soil fungi, and microbes in open and overgrown wetlands. We also tested the effect of encroachment, eutrophication, and soil moisture on indicators of conservation value (red‐listed species, indicator species, and uniqueness). Plant and bryophyte species richness peaked at low soil fertility, whereas soil fertility promoted soil microbes. Soil fungi responded negatively to increasing soil moisture. Lidar‐derived variables reflecting the degree of scrub cover had predominantly positive effects on species richness measures. Conservation value indicators had a negative relationship to soil fertility and a positive to encroachment. For plant indicator species, the negative effect of high nutrient levels was offset by encroachment, supporting our hypothesis of competitive release under shade. The positive effect of soil moisture on indicator species was strong in open habitats only. Nutrient‐poor mires and meadows host many rare species and require conservation management by grazing and natural hydrology. On former agricultural lands, where restoration of infertile conditions is unfeasible, we recommend rewilding with opportunities for encroachment toward semi‐open willow scrub and swamp forest, with the prospect of high species richness in bryophytes, fungi, and soil microbes and competitive release in the herb layer.


| INTRODUC TI ON
Open fens and meadows are characteristic wetland habitats listed on the EU Habitats Directive and targets for conservation (Council Directive 92/43/EEC, 1992). They are species-rich and host large numbers of rare and threatened species (Bedford & Godwin, 2003;Grootjans et al., 2006;van Diggelen et al., 2006;Wassen et al., 2005). Since the mid-20th century, 80% of European wetlands have been degraded or lost due to, for example, encroachment following abandonment of traditional extensive grazing, and eutrophication (Joyce, 2014;Middleton et al., 2006;Verhoeven, 2014).
Scrub encroachment is part of the natural succession process; open habitats grow into a late successional forest in the absence of disturbances, such as lightning-ignited fire, flooding, and grazing (e.g., Bond et al., 2005;Van Wieren, 1995;White, 1979). However, because of human interference, natural disturbances have diminished overall (e.g., Brunbjerg et al., 2014;Middleton et al., 2006;Scholes & Archer, 1997) and the resulting succession has caused widespread scrub encroachment across habitat types and biomes from savannas and steppes to arctic tundra (Naito & Cairns, 2011). In Europe, an increase in vegetation density in the period 2001-2015 has been documented and the vegetation change is likely caused by woody regrowth after the abandonment of livestock grazing (Buitenwerf et al., 2018). Likewise, in Denmark, 17% of the area registered as a meadow in 1992 has now undergone encroachment (Levin & Nainggolan, 2016), however, most of the historical encroachment has happened in the period 1945(Finderup Nielsen et al., 2021.
The pattern is likely to be the same for fens and meadows. When fens and meadows are overgrown with scrub, they lose their legal EU Habitats Directive protection until the scrub eventually grows into late successional swamp forest, which is also protected by the directive (bog woodland 91D0 or Alluvial forests with Alnus glutinosa and Fraxinus excelsior 91E0). In Denmark, nearly 27 million € are spent annually on agri-environmental supplements for livestock grazing and mowing in nature areas to combat encroachment and conserve open habitats (Ministry of Food and Environment, 2015). Despite the effort, only approximately 20% of the semi-natural grasslands, including wet meadows and moors, are currently under active management . This means that most mires are abandoned and subject to free succession. Besides the abandonment of extant fens and mires, many historical fens and meadows have been actively drained, fertilized, and ploughed and are today arable fields and leys. As part of the green transition, a large share of this lowlying farmland is projected to be abandoned and rewetted to avoid further carbon loss from the organic soils. While abandonment from agriculture implies a potential for biodiversity, these areas often have large nutrient pools and strongly modified hydrology due to decades of agricultural use. Eutrophication is a threat to speciesrich open meadow plant communities due to severe asymmetric competition among plant species increasing with high soil fertility (Grime, 1973;Wassen et al., 2005). However, increased shading from encroachment may be hypothesized to relax competition for light among herbs and reduce the competitive exclusion in the field layer as compared to open meadows. Grazing may also partly counterbalance the negative effects of eutrophication (Brunbjerg et al., 2014), but is unlikely to fully compensate (Ejrnaes et al., 2006). Moreover, the full positive effect of disturbances may depend on the restoration of natural hydrology (Kołos & Banaszuk, 2013. The combined effects of nutrients, hydrology, and disturbance regimes in restored wetlands are difficult to predict, but recent studies indicate that the transformation from arable fields to wetlands often fails to restore the species-rich vegetation consisting of stresstolerant forbs and bryophytes characteristic for wetlands (Baumane et al., 2021;Kreyling et al., 2021).
Wetland-restoration success is often evaluated on the basis of plants and birds, while important knowledge obtained from other organism groups, for example, arthropods and fungi, is ignored. In fact, the diversity of heterotrophic organisms, such as arthropods and fungi, is expected to increase with the structural complexity of vegetation and diversity of carbon sources in ecosystems Elton, 1966;Pihlgren & Lennartsson, 2008). In a recent large-scale study, the presence of a shrub layer was the most important variable explaining variation in species richness of fungi and arthropods (Brunbjerg et al., 2020). Heterotrophic organisms gain from the increased biomass following encroachment, as shrubs provide resources and habitats for a large suite of species including herbivores, pollinators, decomposers, and epiphytes (Bruun et al., 2022).
In meadows, eutrophication and groundwater level may control the encroachment process and yield different vegetation structures at different combinations of hydrology and eutrophication. Field observations had led us to hypothesize that thickets and woodlands on wet, nutrient-poor soils grow more heterogeneous in structure, leaving many canopy openings, compared to thickets and woodlands on more nutrient-rich and/or less wet soils. This complexity of vegetation structure can be measured using light detection and ranging (lidar), which is a cost-effective way of gaining fine-resolution data on vegetation structure as compared to field measurements (Lefsky et al., 2002) and which has been shown to capture aspects of vegetation structure that are important and otherwise overlooked for biodiversity (Moeslund, Zlinszky, et al., 2019;Thers et al., 2019). A range of variables representing vegetation structure can be derived from lidar, although the translation to and correlation with well-understood properties is not always straightforward. Despite this, lidar has been used to get detailed information on shrub biomass and cover and has been found to be a "promising tool to map and monitor grassland shrub dynamics at the landscape scale" (Madsen et al., 2020).
In this paper, we investigate the variation in biodiversity along gradients of soil moisture, soil fertility, and scrub cover. We further test the hypothesis that the occurrence of indicators of high conservation value can be promoted by allowing encroachment to take place-especially in restored wetlands on highly eutrophic former agricultural soils. We suggest two mechanisms for such a positive effect: (a) shrubs and trees provide habitat and food resources for large numbers of heterotrophic species, (b) a shrub and tree layer may invoke a competitive release in the herb layer reducing the competitive exclusion of typical wetland plants.

| ME THODS
As part of the present study, we conducted field inventories at 44 wetland sites. The data collection was designed to supplement data collected in the Biowide project, a nation-wide survey of biodiversity in Denmark . Biowide included a total of 130 study sites (40 × 40 m), of which we included all 58 sites evaluated as moist or wet based on plant species composition and soil moisture measurements. The Biowide sites varied in woody species cover from open vegetation over heterogeneous scrub to closedcanopy forest. The Biowide sites also varied in nutrient status from infertile to fertile soils, but were selected to foremost include natural and semi-natural habitats. The additional 44 sites were chosen to increase data coverage of former agricultural land, semi-natural meadows, and agriculturally improved meadows, as well as different levels of scrub encroachment. We did a stratified selection of sites according to succession/light availability (open, tall herb, shrubs, closed canopy), nutrient status (high/fertilized, mid, low/ natural), and soil moisture (moist, wet). Sites were located across Denmark, preferably with a minimum distance of 500 m between sites (except two set of sites, where distances were 252 m and 491 m, Figure 1a). The geographic dispersal of sites ensured sites from both carbon-rich soils, clay, and sand. Each site (40 × 40 m) consisted of four 20 × 20 m quadrants, each with a 5 m circular plot in the center (Figure 1b).
To illustrate the coverage of the soil moisture, nutrient, and encroachment gradients covered by the combined data, we compared site mean Ellenberg F, N, and L values (Ellenberg et al., 1992) for all 5 m circle plots with a reference dataset from national monitoring, using identical 5 m circular plots (59,227 sites from agricultural, semi-natural, natural open, and forest vegetation http://www.natur data.dk) (Nygaard et al., 2017). Mean Ellenberg values were calculated for plots with more than five species present. In scatterplots of plot-mean Ellenberg values, 95 percentile convex hull polygons were drawn for the reference dataset as well as the Biowide and wetland dataset to visually evaluate the representativity of our data (Appendix A).

| Biowide data collection
2014 with a few early spring vascular plant species added in May 2015 . We aggregate the four 5 m circle plots and additional species list to obtain a site species list for analyses.

| Collection of soil eDNA data
As alternative measures of biodiversity, we used the richness of operational taxonomic units, that is, OTUs (Blaxter et al., 2005) of fungi and (eukaryotic) soil microbes obtained from metabarcoding of soil-extracted eDNA (Frøslev et al., 2017

| Wetland sites data collection
All additional data collection specifically for the present study was done according to Biowide protocols , with the exception of the following: (1)

| DNA extraction and metabarcoding
For Biowide and wetland soil samples, DNA was extracted and subjected to eDNA metabarcoding through DNA extraction, PCR amplification of genetic marker regions (DNA barcoding regions), and massive parallel sequencing on the Illumina MiSeq platform as described in Brunbjerg et al. (2019). For this study, we used highthroughput sequencing data from marker genes amplified with primers targeting eukaryotes (mainly soil microbes) (18S) and fungi (ITS2).
OTU tables were constructed following the overall pipeline in Frøslev

| Lidar-based measures
We calculated 23 lidar variables to represent encroachment in the 102 sites. We used the same procedure as in Valdez et al. (2021).

| OTU richness
As alternative measures of biodiversity, we used the richness of operational taxonomic units, that is, OTUs (Blaxter et al., 2005) of fungi and soil microbes from metabarcoding of soil eDNA (Frøslev et al., 2017. Classical data collection of fungi is time consuming and OTU richness has been found to resemble classical observed species richness at least for groups of macrofungi that are feasible to include in field inventories . We used OTU richness of fungi and soil microbes as response variables to reflect diversity of species groups not monitored otherwise in this project.  Commission, 2007). Common to these indicator species is a preference for infertile habitats (low Ellenberg N and high Grime's S values, Grime, 1979).
Biotic uniqueness: Uniquity is a scale-dependent metric of biodiversity reflecting how unique the biodiversity of a given site is compared to the gamma diversity across the containing region or collection of sampled sites . Uniquity can be calculated based on both observational data as well as non-annotated DNA data (e.g., OTUs) and hence can reflect both species uniqueness and genetic uniqueness. Contrary to other biodiversity metrics, uniquity accounts for sampling bias and spatial scale. Due to the built-in weighting method, uniqueness of non-annotated DNA-data can be calculated corresponding to, for example, red-listed species . Here, we calculated fungi and soil microbe uniquity according to Ejrnaes et al. (2018) in order to reflect the unique site con-  Ellenberg et al., 1992) has been found to reflect eutrophication in wetlands and be highly correlated with the number of typical species in fens (Andersen et al., 2013). For each site, we calculated mean Ellenberg N and Ellenberg R values (plant-based bioindication of nutrient status and soil pH, respectively) (Ellenberg et al., 1992). The Ellenberg nutrient ratio was used to reflect eutrophication and the idea of the ratio is to account for the fact that natural nutrient availability in wetlands increases with pH. To avoid circularity in analyses, plant species included in the plant-based con- Standardized plant species richness was used as co-variable in models using OTU richness as response. Negative binomial errors were used if overdispersion was detected in Poisson models (Hilbe, 2011). We allowed for interaction between lidar variables and nutrient ratio, and lidar variables and soil moisture, respectively. We included a quadratic term of nutrient ratio and moisture variables if the full model significantly improved according to the ΔAIC < 2 criterion (Burnham & Anderson, 2002). The residuals of full models were checked for model misfit and overdispersion and spatial autocorrelation using correlograms from the R package ncf (Bjørnstad, 2020). Because of the large number of explanatory variables we used stepwise forward selection using the ΔAIC < 2 criterion (Burnham & Anderson, 2002).
We used the variation inflation factor (VIF) to test for co-variability among selected explanatory variables and accepted values <3.

| Conservation models
In order to test the effect of encroachment on conservation value,  (Hilbe, 2011). We allowed for interaction between encroachment and nutrient ratio, and encroachment and soil moisture, respectively. We included a quadratic term of the explanatory variables if the full model significantly improved according to the ΔAIC < 2 criterion. The residuals of full models were checked for model misfit and overdispersion and spatial autocorrelation using correlograms from R package ncf (Bjørnstad, 2020). We used backwards elimination of explanatory variables using the ΔAIC < 2 criterion (Burnham & Anderson, 2002) to reduce full models to final models.
We checked if soil eDNA data for the two datasets could be

| RE SULTS
The study sites covered the full gradient in nutrient availability of the reference data, but only the wetter part of the reference moisture gradient as expected (Appendix A). The encroachment gradient (Ellenberg L) was also almost fully covered. The wetland dataset sup-

| Richness models
We found a negative effect of nutrient ratio on plant species richness. However, only c. 7% of the variation in plant species richness could be accounted for by the model ( Table 1)

| Testing the effect of encroachment on conservation interest
The explanatory strength of conservation models ranged from 9.5% for the soil fungal uniquity model to 41.3% for the model for indicator species (Table 2). Across all response variables, high degree Encroachment is often considered a threat to open-landscape biodiversity (Ratajczak et al., 2012;Stoate et al., 2009). In a conservation management perspective, focus is often on maintaining early successional vegetation, for example, by grazing and mowing of fens (van Diggelen et al., 2015) to ensure favorable conditions for especially rare plant species sensitive to encroachment (Bart, 2021).
Traditionally, there has been a focus on vascular plants, when defining habitat types, evaluating conservation status, and planning management, for example, within the framework of the EU Habitats Directive . However, this plant-focus may create a biased perception of the effect of encroachment on biodiversity as this effect at least depends on the habitat type and the response group in question (Eldridge et al., 2011). To our surprise, we did not find a negative effect of encroachment in our models for red-listed vascular plant richness or indicator plant richness. This is not to say that light-demanding species are not replaced by invading shrubs, but either the loss is relatively weak or losses are off-set by gains of equally rare species. For the indicator plant species, the highest richness is observed in infertile and open habitats (Figure 3) pointing to the need for protecting these against encroachment.
We found a positive effect of encroachment on bryophyte species richness and presence of red-listed bryophytes. Bryophyte richness was higher in sites with relatively high vegetation cover and low lidar amplitude. A high vegetation cover (measured by lidar with 5 points/m 2 ) does not necessarily mean that the vegetation is dense (recall that the lidar was recorded during leaves-off and hence it is unlikely to capture the herb layer reliably). Instead, it likely means that bryophyte richness is highest where the vegetation includes TA B L E 1 Modeling output for GLM negative binomial models using vascular plant species richness, bryophyte species richness, fungal OTU richness, and soil microbe OTU richness as response variables and nutrient ratio (Ellenberg N/R), soil moisture, and the set of 23 lidar variables as explanatory variables.  small shrubs and trees, which also seem to be the case for the richness of red-listed bryophytes in our study.
Heterotrophic organisms are likely to benefit from the expansion of niches linked to build-up and diversification of organic carbon following encroachment with shrubs and trees (Brändle & Brandl, 2001;Bruun et al., 2022). that in nutrient-poor wetlands, trees contribute to ecospace expansion , maybe in the form of substrate (i.e., falling leaves or root sap). Lastly, soil microbe richness was higher at TA B L E 2 Modeling output for models using number of red-listed vascular plants (logistic), red-listed bryophytes (logistic), number of indicator species (GLM negative binomial), fungal uniquity (Gaussian), and soil microbe uniquity (Gaussian) as response variables and nutrient ratio (Ellenberg N/R), soil moisture, and encroachment as explanatory variables.  For soil fungi and microbe communities, we found more unique assemblages with encroachment, but for indicator species of vascular plants in contrast, we found the classic peak of high species richness in open, nutrient-poor fens (Wassen et al., 2005). Others have also found complex richness responses to shrub encroachment, including a hump-shaped relationship (Kesting et al., 2015), indicating that encroachment of scattered shrubs in open grasslands may cause increased habitat heterogeneity which benefit species richness, while complete overgrowth will lead to reduced vascular plant species richness on the scale of small sample plots (Dierschke, 2006;Galvánek & Lepš, 2008;Ratajczak et al., 2012;Teleki et al., 2020), probably due to light extinction and leaf litter cover inhibiting seedling establishment (Jensen & Schrautzer, 1999).

Red-listed plants
Eutrophication is well-documented as a major threat to freshwater meadow biodiversity. Eutrophication causes a shift in species composition from slow-growing, light demanding vascular plants and bryophytes to more competitive and fast-growing species (Bobbink, 2004;Bobbink et al., 1998;Hogg et al., 1995)-a more rapid shift than the vegetation changes due to natural succession (Hogg et al., 1995). Our results are aligned with the negative effect of soil fertility for both species richness of vascular plants and bryophytes, but show a more complex interaction for soil microbial OTU richness. Hence, the diversity of soil microbes increased with eutrophication in open meadow sites, possibly reflecting increased available carbon to use as substrate. This effect is absent from encroached sites, possibly because the litter and rhizosphere of shrubs and trees add diverse carbon sources, irrespective of soil fertility.
Species of conservation interest and uniqueness of soil fungi similarly showed a negative response to soil fertility. The number of redlisted plant and bryophyte species decreased with increasing soil fertility. Negative effects of eutrophication has been found to be more severe for rare species due to their initial low abundance-at least in grasslands and wetlands (Clark & Tilman, 2008).
We hypothesized that in eutrophic sites, competitive release (Keddy & Maclellan, 1990) may be a positive consequence of scrub encroachment and the resulting vertical differentiation of vegetation layers. The competitive release hypothesis is underpinned by the notorious depauperate plant species richness in eutrophic herbaceous vegetation due to asymmetric competition for light and nutrients (e.g., Crawley et al., 2005). We found the mentioned strong negative effect of eutrophication on plant species richness in open herb communities, but a much smaller effect under canopy cover for the number of indicator species (Figure 3), supporting the hypothesis of competitive release.
Rewetting and recreating natural hydrology is a well-established management recommendation for fen and meadow systems (Kołos & Banaszuk, 2013. However, we did not find a general positive effect of soil moisture on species richness and indicators of conservation value in our study, the effect was only positive for indicator species in open habitats but not after encroachment. We suspect that the effect of soil moisture can be partly confounded with both eutrophication and encroachment because leached nutrients in the watershed are transported with the water and released into the wetland communities and also the wettest areas tend to be abandoned first and generally exhibit heavier encroachment than less wet sites.
Anecdotic evidence from our own surveys of aerial photographs indicates that willow scrub and alder swamps predominantly occur in the wettest parts of river valleys, for example, in places where historical small-scale peat extraction has left inundated pits and rendered the tract unsuitable for cultivation.
In our study, plant species richness and the number of red-listed plants were solely affected by soil fertility and the models only had limited predictive power (c. 7% and c. 14%, respectively). The predictive power of these models may seem low when compared to c.
60% explained variation in plant species richness in the full Biowide dataset in Brunbjerg et al. (2020). Although the two studies are not directly comparable, several explanations may be suggested for this discrepancy, for example, our encroachment variable is rather crude and also variables representing disturbance (e.g., grazing) and historical events related to former cultivation not included in this study could possibly be important for the plant species richness in our wetland sites.
Because of the need for reducing greenhouse gas emissions (the Paris agreement, United Nations, 2015), for example, by agricultural abandonment of organic soils/peat lands to decrease CO 2 -emissions, Denmark is now planning to abandon 100.000 ha of carbon-rich lowlying cultivated areas. These areas are likely to have high nutrient status from former cultivation and restoring natural low nutrient status as a basis for developing protected, open habitat types are costly and tedious. Instead, allowing encroachment when abandoning is expected to further increase CO 2 -sequestration because of the accumulation of carbon in shrub biomass and will furthermore permit synergies between climate and biodiversity goals. While setting aside cultivation of these potential wetlands implies a great potential for ecological restoration, our study shows that notably, earlier eutrophication caused by decades of arable farming will almost inevitably hamper the restoration target of species-rich meadows and fens.
Based on our results, we recommend to combine a relaxed attitude to encroachment with reintroduction of natural disturbances (e.g., widespread rewilding of large herbivores) in order to promote semi-open scrub and woodland communities. Scattered bushes and thickets are natural elements in grazing systems, as many shrub species are vigorous resprouters, for example, Salix species (Klimkowska et al., 2010).
In areas not suitable for year-round grazing, so-called "passive rewilding" (i.e., natural processes are allowed to restore themselves, Svenning et al., 2016)

ACK N OWLED G M ENTS
This work is a contribution to the project "Mapping, restoration and management of groundwater-fed fens and springs" funded by Aage V. Jensen Naturfond. The Biowide project was supported by a grant from the Villum Foundation (VKR-023343).

A PPEN D I X B
Pictures of nine of the 102 sites and their approximate position in an ecological space of encroachment represented by vegetation density at 3-10 m (x-axis) and productivity/nutrient status represented by Ellenberg N (y-axis).

A PPEN D I X C D et ail e d info r m at io n o n m o l e cu l a r d at a DNA extraction (44 study sites)
DNA was extracted from 4 g of soil using the PowerMax Soil DNA Isolation kit (Qiagen), following the manufacturer's procedure.
DNA extract was purified using the PowerClean DNA Clean-Up Kit (Qiagen), and DNA was normalized to 1 ng/μl after fluorometric quantification using the Qubit™ dsDNA HS Assay Kit (Thermo Fischer). The 44 samples were extracted in smaller batches with one negative control for each batch.

Sequence processing
The bioinformatic processing of the sequence data followed the strategy outlined in Brunbjerg et al. (2019). Demultiplexing of samples was done with custom scripts that keeps R1 and R2 separate for DADA2 processing , and is based on Cutadapt (Martin, 2011)-and also Sickle (Joshi & Fass, 2011) for fungi-as in Frøslev et al. (2019) searching for a sequence pattern matching the full-length combined tag and primer allowing for no errors, and removing possible remnants of the other primer at the 3′ end. We used DADA2 v 1.8 (Callahan et al., 2016) to identify OTUs (also known as exact amplicon sequence variants, ESVs) and for removal of chimeras (bismeras). For eukaryotes, the OTU tables were used directly.
The fungal data were clustered at 98.5% and filtered to contain only ingroup data-that is, kingdom Fungi. For the fungal dataset, taxonomy was assigned by matching against the UNITE database (Nilsson et al., 2019), and for the eukaryotic data using a custom script based on BLASTn searches against genbank (Altschul et al., 1990).
Sequence data and analytical documentation can be obtained by contacting the first author.

A PPEN D I X D
Lidar-based explanatory variables and interpretation. If the standard deviation of a variable was calculated, in addition to its mean, the variable is denoted with an asterisk.

Explanatory variable Explanation
Light/heat Potential solar radiation* Potential solar radiation input only considering the terrain Adjusted potential solar radiation* Same as above but adjusted for vegetation cover Amplitude Succession state (high amplitude = high reflectance/ flat and simple vegetation, low amplitude = low reflectance/complex vegetation). The measure is corrected for aircraft type and seasonality, (see Valdez et al., 2021) Echo ratio* Canopy cover and complexity

A PPEN D I X E
Indicator species considered moderately to very sensitive toward habitat changes as defined by Fredshavn et al., 2010. The list of indicator species was developed to indicate favorable conservation status cf. the Habitats Directive (European Commission, 2007). Species name and family are given.

A PPEN D I X F
Variation in ecological conditions and species richness across the 102 sites (mean and standard deviation are given). Fine-scale terrain roughness 0.01 0.00